Seismic sparse-layer reflectivity inversion using basis pursuit decomposition

Basis pursuit inversion of seismic reflection data for reflection coefficients is introduced as an alternative method of incorporating a priori information in a seismic inversion process. The inversion is accomplished by building a dictionary of functions, representing seismic reflection responses, and constituting the seismic trace as a superposition of these responses. Basis pursuit decomposition finds a sparse number of reflection responses that sum to form the seismic trace. When the dictionary of functions is chosen to be a wedge model of reflection coefficient pairs convolved with the seismic wavelet, the resulting reflectivity inversion is a sparse-layer inversion, rather than a sparse-spike inversion. Synthetic tests show that sparse-layer inversion using basis pursuit can better resolve thin beds than a comparable sparse-spike inversion can. Application to field data indicates that sparse-layer inversion results in potentially improved detectability and resolution of some thin layers and reveals apparent stratigraphic features that are not readily seen on conventional seismic sections. INTRODUCTION In conventional seismic deconvolution, the seismogram is convolved with a wavelet inverse filter to yield band limited reflectivity. The output reflectivity is bandlimited to the original frequency band of the data so as to avoid blowing up noise at frequencies with little or no signal. It has long been established (e.g., Riel and Berkhout; 1985) that sparse seismic inversion methods can produce output reflectivity solutions that contain frequencies that are not contained in the original signal without necessarily magnifying noise at those frequencies. It is well known (e.g., Tarantola, 1987) that applying valid constraints in seismic inversion can stabily increase the bandwidth of the solution. However, incorporation of a priori information in the reflectivity inversion of seismic traces can be problematical. A common way of incorporating prior knowledge is to build a starting model biased by that information and letting the inversion process perturb the initial starting model and converge to a solution (e.g., Cooke and Schneider; 1983). The individual layers represented in the starting model can have hard or soft constraints assigned. This kind of method can work very well when the starting model is close to the correct solution. Typically, the starting model is obtained by spatially interpolating well logs along selected horizons. Unfortunately, these horizons must be picked on the original seismic data. If waveform interference patterns change laterally, horizon picks on a constant portion of a waveform (typically chosen to be peaks, troughs, or zero crossings) can be in error, resulting in an incorrect starting model and possibly an erroneous inversion. Similarly, if velocities and/or impedances for the inversion interval change laterally in a manner different from that resulting from the interpolation procedure, interpolated well logs may again be significantly in error and the inverse process may converge to the wrong minimum. It may be possible to ameliorate these problems to some degree using a Monte Carlo approach, but this cannot correct the fundamental non-uniqueness of the process that may cause a number of minima other than the correct one to have similar error. A means of biasing the results towards expected reflectivity patterns is needed without relying on possibly erroneous manual interpretations or spatial interpolations. Nguyen and Castagna (2010) used matching pursuit decomposition (MPD) to decompose a seismic trace into a superposition of reflectivity patterns observed in and derived from existing well control. Matching pursuit decomposition (1) correlates a wavelet dictionary against a seismogram and finds the location, scale (i.e.; center frequency), and amplitude of the best fit wavelet, (2) subtracts the best fit wavelet and records its characteristics in a table, and (3) repeats the processs on the residual trace until the residual energy falls below a selected threshold. For spectral decomposition the spectra of the best-fit wavelets from the table are summed as a function of time to form the time-frequency analysis. For seismic inversion, the wavelet dictionary consists of seismic reflection patterns derived from well logs which are matched to the seismic trace. In effect, pattern recognition is used to recognize seismic patterns derived from well logs resulting in what is equivalent to a data adaptive starting model. This ensures that the starting model is not misaligned with the seismic data to be inverted as may occur when well logs are spatially interpolated. The matching pursuit method has some limitations, especially when the dictionary elements are not orthogonal. Nguyen and Castagna (2010) consequently found lateral instability in the raw MPD solution and had to employ re-projection in order to obtain laterally stable results. Wang (2007, 2010) developed a multichannel MPD (MCMP) spectral decomposition assuming some degree of lateral coherence in order to improve the uniqueness and spatial continuity of matching pursuit spectral decomposition. As applied to the seismic inversion problem, such an approach has potential limitations as it (1) could result in some loss of spatial resolution when the geology is, in fact, not continuous, and (2) can still exhibit hard lateral jumps when the path-dependent MPD algorithm switches from one initial match to another. .For seismic inversion purposes, a more laterally stable approach that operates on a single trace at a time, without the need for a posteriori selection of laterally consistent possible solutions, would be advantageous. Basis Pursuit Decomposition (BPD) has a number of advantages over MPD (Chen et al., 2001); it handles interferences between dictionary elements better, is computationally more efficient and, by introducing a sparsity norm and regularization parameter into the objective functions, it can exhibit good lateral stability even when dictionary elements are not orthogonal. In this study, we investigate the use of BPD to perform seismic inversion. We use the algorithm of Chen et al. (2001) to decompose the seismic trace with a nonorthogonal wavelet dictionary consisting of even and odd thin layer seismic responses. The output reflectivity series is then formed by summing reflection coefficient pairs that are shifted and multiplied by scalars with translations and coefficients output by the basis pursuit algorithm. The basis pursuit algorithm is an L1 norm optimization that was originally developed as a compressive sensing technique. The solution consists of three parts: (1) linear programming transforms the L1 optimization problem into a constrained least square problem; (2) duality theory sets up an array containing primal, dual and gap equations to be solved; and (3) a primal-dual log-barrier algorithm implements a GaussNewton step workflow to solve the equations array. Much of the potential advantage of BPD over MPD is related to the fact that BPD finds a single global solution, whereas MPD is a path-dependent process. As MPD iteratively subtracts matched wavelets from the seismogram, the order in which wavelets are subtracted may vary greatly between similar seismic traces. A change in the order in which the wavelets are subtracted can result in an entirely different solution. Slight variations in adjacent input traces (as could result from noise) could cause geologically unreasonable jumps from one solution to another. Furthermore, for interfering wavelets, the location and center frequency of the best wavelet match found by MPD may not correspond to the time location or center frequency of the reflecting seismic wavelets. If a wavelet of the wrong frequency is subtracted, nearby wavelets found by MPD must compensate for this error; often by being wrong in the opposite direction, and also can result in the addition of small satellite wavelets to the solution that are needed to minimize the mismatch. In data compression applications this results in less than optimal compression, but is not necessarily a serious problem. In spectral decomposition using MPD however, this can cause time asymmetry in the timefrequency decomposition of thin layers, even when the input is time symmetrical. This results in serious lateral instability in spectral decomposition and inversion using traceby-trace MPD (Wang, 2010; Nguyen and Castagna, 2010). Our investigation of the use of basis pursuit inversion (BPI) to invert a seismic trace for a reflectivity series which can be integrated to give the band-limited seismic impedance model, is an attempt to take advantage of the ability of MPD to incorporate a priori information into the wavelet dictionary in an inherently more stable manner. As a first step, we study the use of a very basic reflection pattern dictionary consisting of a wedge model of reflection coefficient pairs convolved with the seismic wavelet. As any pair of reflection coefficients can be represented as the sum of odd and even reflection coefficient pairs (e.g., Puryear and Castagna, 2008), there will be two dictionary elements (odd and even) for each thickness represented in the wedge model and the number of elements in the dictionary will be twice the number of thicknesses. The thickness range built into the dictionary then becomes a constraint on the possible outcomes. The direct use of more complex well log derived reflectivity patterns in the BPI dictionary is left as an objective of future research. The use of BPI with a wedge dictionary is essentially a sparse-layer inversion. It can be argued that such an inversion may have different ability to resolve thin beds, than a sparse-spike type of inversion (e.g., Oldenburg et al.;1983) that imposes some specified degree of sparsity and thus will indirectly place some limitation on the spacing of reflections. In this paper we compare BPI sparse-layer inversion with conventional sparse-spike inversion on synthetic and real data to assess the relative ability of BPI to resolve thin beds and reveal fine stratigraphic features. BASIC THEORY For simplicity, the forward model is assumed to be a simple convolution of a stationary seismic wavelet and the reflectivity. The seismic trace s(t) is thus given by